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Abstract 

Living systems are forced away from thermodynamic equilibrium by exchange of mass and 
energy with their environment. In order to model a biochemical reaction network in a 
non-equilibrium state one requires a mathematical formulation to mimic this forcing. We 
provide a general formulation to force an arbitrary large kinetic model in a manner that 
is still consistent with the existence of a non-equilibrium steady state. We can guarantee 
the existence of a non-equilibrium steady state assuming only two conditions; that every 
reaction is mass balanced and that continuous kinetic reaction rate laws never lead to a 
negative molecule concentration. These conditions can be verified in polynomial time and 
are flexible enough to permit one to force a system away from equilibrium. In an expository 
biochemical example we show how a reversible, mass balanced perpetual reaction, with 
thermodynamically infeasible kinetic parameters, can be used to perpetually force a kinetic 
model of anaerobic glycolysis in a manner consistent with the existence of a steady state. 
Easily testable existence conditions are foundational for efforts to reliably compute non- 
equilibrium steady states in genome-scale biochemical kinetic models. 
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1. Introduction 

There are various approaches for simulation of biochemical network function [19j. In 
principle, an ideal approach would accurately represent known physicochemical principles 
of reaction kinetics, tailored with kinetic parameters specific to a particular organism. 
However, when modeling genome-scale biochemical networks, one's choice of modeling 
approach is also shaped by concerns of computational tractability. One of the main reasons 
that flux balance analysis [36\ \TT\ 131} [2^ [23] has found widespread applications in genome- 
scale modeling is that the underlying algorithm is typically based on linear optimization. 
In general, industrial quality software implementations of linear optimization algorithms 
are guaranteed to find an optimal solution, if one exists, or otherwise give a certificate that 
the problem posed is infeasible. In the process of iterative model development, one may use 
flux balance analysis to test if a stoichiometric model, obtained from a draft reconstruction, 
actually admits a steady state flux (reaction rate) [33j. If not, various algorithms have been 
developed that help to detect missing reactions |33j or stoichiometric inconsistencies [14J 



(discussed further in Section 2.1). 



Flux balance analysis predicts fluxes that satisfy steady state mass conservation but not 
necessarily energy conservation or the second law of thermodynamics [H ^I7\ I13| . Whilst 
the set of steady state mass conserved fluxes includes those that are thermodynamically 
feasible, additional constraints are required in order to guarantee a flux that additionally 
satisfies energy conservation and the second law of thermodynamics [12] . Within the set of 
steady state fluxes that satisfy mass conservation, energy conservation and the second law 
of thermodynamics, there is a subset that additionally satisfy various reaction kinetic rate 
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laws [8] . The satisfaction of such kinetic constraints is important for accurate representation 
of various biochemical processes where the abundance of a molecule affects the rate of a 
reaction, e.g., allosteric regulation [161 [18]. 

There are various algorithmic barriers to genome-scale kinetic modeling that preclude 
satisfaction of all of the aforementioned thermodynamic and kinetic constraints, without 
resorting to rate law approximation. Apart from the open algorithmic challenge to de- 
velop an algorithm for efficient computation of thermodynamically and kinetically feasible 
steady states in genome-scale kinetic models, there is also the challenge of mathematically 
expressing the necessary and sufficient conditions for existence of such a steady state in a 
manner that can be efficiently tested given a genome-scale model. 

The quantitative study of chemical reaction kinetics has a long history, beginning per- 
haps with Ludwig Wilhemy's 1850 discovery that the rate of a chemical reaction is propor- 
tional to the concentrations of consumed substrates [39] . This fundamental law of chemical 
kinetics is well known to chemists and biochemists alike. However, due to imperfect in- 
heritance of knowledge, there are other useful facts, established generations ago, that are 
slipping from the consciousness of chemists and biochemists alike - as measured by con- 
temporary citations. Even for a multifaceted paper with many citations, those citations 
can be due to a historically facet, not necessarily the most useful facet in a contemporary 
sense. A case in point is the 1962 paper by James Wei entitled "Axiomatic treatment of 
chemical reaction systems" [37]. This excellent paper has been cited 65 times, but only 
twice in the last 10 years and infrequently by theoretical biochemists. The vast majority 
of citations refer to Wei's treatment of the stability of chemical reaction systems with Lya- 
punov functions. Exceptionally, the importance of Wei's result, concerning the conditions 
for existence of non-equlibrium steady states, is realized, e.g., in 1976 Perelson used Wei's 
result to illustrate the danger inherent in concluding results from mathematical models of 
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systems of chemical reactions that do not conserve mass [25]. 

Herein we build upon Wei's work that establishes sufficient conditions for the existence 
of at least one steady state reaction flux and molecule concentration for a broad class 
of chemical kinetic models. This class of kinetic model includes all those networks with 
exclusively mass balanced chemical reactions, where the kinetics rate laws are such that 
molecule concentration can never be a negative quantity. This class of reaction network 
includes all biochemical reaction networks. The conditions for existence may be easily 
tested by a trivial check on the kinetic rate law formulation for each reaction, together 
with a test of stoichiometric consistency using linear optimization |14j . 

In Section [2] we introduce some mathematical definitions of pertinent chemical reac- 
tion network concepts. Section [3] states and proves theorem concerning the existence of 
a non-negative steady state molecule concentration (vector) for mass conserved elemen- 
tary reaction kinetics. This theorem and proof follow the more general case outlined in 
broad strokes by Wei j^. Based on citation history, Wei's result seems to be lost upon 
the biochemical modeling community. In Section [4j we illustrate for the first time, the 
utility of Wei's existence theorem for modeling non-equlibrium steady states with various 
examples, including a kinetic model of anaerobic glycolysis in Trypanosoma brucei. Finally 
we summarise and attempt to place this work in the context of established mathematical 
approaches to model biochemical reaction networks. 

2. Chemical reaction networks 
2.1. Stoichiometry 

Consider a biochemical network with m molecules and n elementary reactions. An 
elementary reaction is one for which no reaction intermediates have been detected or need to 
be postulated in order to describe the chemical reaction on a molecular scale. It follows that 
the reaction stoichiometry is sufficient to define the molecularity of the molecules involved in 
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the reaction. One may combine elementary reactions together to form a composite reaction. 
One can define the topology of the resulting hypergraph using a generalized incidence 
matrix, S G Z™"'", where S is always singular and typically r = rank(S) < m < n for 
large biochemical networks. Each row in this stoichiometric matrix represents a particular 
molecule, e.g., glucose, whilst each column represents a reversible biochemical reaction. 
We assume that all biochemical reactions are indeed reversible [20j. For each reversible 
reaction, convention dictates one direction be designated forward and the other reverse. 
With respect to the forward direction, for all i = 1 . . . m and j = 1 . . . n, Sij < if molecule 
i is a substrate in a reaction, meaning that it is consumed by the reaction j, Sij > if 
molecule z is a product, meaning that it is produced by a reaction, and Sij = otherwise. 
Typically stoichiometric coefficients are integers reflecting the whole number molecularity 
for a molecule consumed or produced in a reaction. 

Each column of a stoichiometric matrix contains at least one negative coefficient and 
one positive coefficient, reflecting either the chemical conversion of one molecule to another, 
or in multi-compartmental models, the transport of a molecule from one compartment to 
another, i.e., a transport reaction may consume one molecule in a reactant compartment 
and produces one molecule in a different product compartment, even if the molecule is 
physically identical. We assume that each column of S corresponds to one mass conserving 
chemical reaction. A necessary, but insufficient condition for mass balancing is that each 
column of S must have at least one positive coefficient and at least one negative coefficient. 
We say that a chemical reaction is linear when the corresponding column of S contains 
two nonzero coefficients, {—1, 1}. We say that a chemical reaction is bilinear when the cor- 
responding column of S contains three non-negative coefficients, {—1, 1, 1} or {—1, —1, 1}. 
There may be more than one negative (positive) coefficient in a column when a reaction 
involves more than one substrate (product). In reaction networks with composite reactions, 
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nonzero stoichiometric coefficients are typically not of magnitude one. However, even the 
most complicated composite reaction can be decomposed into linear and bilinear reactions. 

Each row of S contains at least one positive coefficient and at least one negative co- 
efficient, reflecting the requirement for at least one reaction to produce and at least one 
reaction to consume each molecule. A stoichiometric matrix for a chemical reaction net- 
work is said to be consistent if each molecule can be assigned a single positive molecular 
mass, without violating mass conservation, and inconsistent otherwise |14tll0j. Mathemat- 
ically, this translates to the existence of at least one strictly positive vector, m S IK>0) in 
the left nullspace of a consistent stoichiometric matrix S"^ • m = 0. Strictly, we could say 
that each row of S corresponds to an isotopically distinct molecular entity in order that 
the corresponding molecular mass be precisely defined, as two otherwise identical molecules 
can have different molecular mass depending on their isotopic label e.g. ^^C vs ^^C glucose. 
However, we shall assume that reaction kinetic parameters are isotopomer invariant so we 
need not be so strict. 

2.2. Reaction kinetics 

Perhaps the simplest reaction kinetic assumption is that a unidirectional reaction rate 
is proportional to the product of the concentrations of each substrate consumed [39] • Let 
us define forward and reverse stoichiometric matrices, F, R G M^q" respectively, where 
Fij denotes the stoichiometry of substrate i in forward reaction j and Rjj denotes the 
stoichiometry of substrate i in reverse reaction j. It follows that the stoichiometric matrix 
is defined by S = — F+R. It is possible for the same molecule to appear as both a substrate 
and a product in the same unidirectional reaction, e.g., an auto- catalytic reaction, so it 
is natural to define S in terms of F and R, rather than the other way around. We we 
may now express elementary kinetics for forward and reverse reaction rates, respectively 
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vj, Vr G M", as 

v/(k/,x) = diag(k/) • exp(F'^ • ln(x)), 

Vr(kr,x) = diag(kr) • exp(R-^ • ln(x)), 
where we assume non-negative elementary kinetic "parameters kjjk^ G IK>0) ^-ncl the expo- 
nential or natural logarithm of a vector is meant component- wise We say that a pair 
of forward and reverse elementary kinetic parameters are thermodynamically feasible when 
they satisfy 

where u° £ M™ is a vector of standard chemical potentials, R is the gas constant and T 
is temperature. Thermodynamically feasible kinetic parameters for all reactions implies 
detailed balance at thermodynamic equilibrium, i.e., vj = |35j . 

We refer to ([T]) as mass action kinetics [28\ only after we have stated our assumption that 
Q also holds for each mass conserving reversible reaction. In the words of Horn & Jackson 
|17j "a kinetic description of chemical reactions in closed systems with ideal mixtures, 
completely consistent with the requirements of stoichiometry and thermodynamics, may 
be obtained by satisfying the following four conditions": 

(a) The rate function of each elementary reaction is of the mass action form. 

(b) The stoichiometric coefficients are such that mass is conserved in each elementary 
reaction. 

(c) The kinetic constants in the rate functions are constrained in such a way that the 
principle of detailed balancing is satisfied. 



^Strictly, it is not proper to take the logarithm of a unit that has physical dimensions. This difficulty 
can be avoided by considering x as a vector of mole fractions rather than concentrations (Eq. f9.93 in [Q). 
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(d) The stoichiometric coefficients are non-negative integers. 

Condition (a) is represented by ([T]) and conditon (b) is satisfied by strict adherence to 
elemental balancing for each reaction during network reconstruction \33\ I34j . Horn Sz 
Jackson [17J considered general mass action kinetics when only (a) was assumed to hold. 
We shall consider mass conserved elementary kinetics where (a) and (b) are assumed to 
hold but (c) and (d) are allowed to be relaxed for a subset of reactions. We shall return to 
this point in the discussion. We take the dynamical equation for mass conserved elementary 
kinetics to be 

r/x 

X = — = S • (K/ • exp(F'^ • In (x)) - • exp(R^ • In (x))) (3) 

where t denotes time, all reactions conserve mass, all reactions are reversible, Kj = 
diag(kj), Kr = diag(kr) and kj,kr € K™q. 

2.3. Concentration non-negativity 

Physically, one would expect that molecular concentration be a non-negative quantity. 
Starting from an initial non-negative concentration, xq > 0, it has been proven math- 
ematically that all subsequent concentrations are non-negative when the evolution of a 
system is subject to elementary reaction kinetics [SljT]. Elementary kinetics is essentially 
non-negative if, for all i = 1 . . . m, > for all Xj G K>o, where Xi and Xi denote the 
j^th component of x and x, respectively. This mathematical formulation can be easily be 
chemically interpreted. Suppose the concentration of molecule A is zero; then irrespective 
of what the non-negative concentration is for all other molecules, the rate of change in 
concentration for molecule A is non-negative. To understand why, recall that in elemen- 
tary kinetics the rate of a unidirectional reaction is always the product of a non-negative 
elementary kinetic parameter and the concentration(s) of the substrate(s), each to the 
power of the absolute value of the corresponding stoichiometric coefficient. If a molecule's 
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concentration is zero then all reactions consuming that molecule have a rate of zero; hence, 
there can be no consumption of that molecule and its rate of change in concentration is 
non- negative. 

Consider the chemical reaction network with the single reaction 

A^B + C. 

Let the forward and reverse reaction rates be given by elementary kinetics Vf = kfU and 
Vr = krbc with kf, kr G M>o, where lowercase refers to concentration of the corresponding 
uppercase molecule. If the initial concentrations of all molecules are non-negative it is 
impossible for any molecule's concentration to go negative because a = implies Vf = 
so no consumption of A would occur and therefore a cannot be negative. Observe that 
the conditions for concentration non-negativity place no constraints on the ratio of forward 
over reverse kinetic parameter, but only that all kinetic parameters be non-negative. 

A further technical restriction is that x be given by a locally Lipschitz continuous 
function, but this is easily satisfied by deterministic formulations of kinetics where the cor- 
responding differential equations are continuously differentiable (Lemma 1, [7|). Starting 
from an initial non-negative concentration, equation ^ constrains all subsequent concen- 
trations to be non-negative as we assume elementary reaction kinetics with non-negative 
kinetic parameters. Even for non- negative but thermodynamically infeasible kinetic param- 
eters, which violate ([2]), it is true that if one begins with non- negative initial concentrations 
then all subsequent concentrations remain non-negative. 

3. Existence of a steady state for mass conserved elementary kinetics 

The following theorem establishes sufficient conditions for existence of a steady state 
concentration for a dynamical system governed by mass conserved elementary kinetics. 
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Theorem 1. Let the dynamical equation for mass conserved elementary kinetics be 



dx 

~dt 



S • (K/ • exp(F^ • In (x)) - • exp(R'^ • In (x))) 



(4) 



where x = x(t) G is molecule concentration at time t > 0, x G M™ is the time 
derivative of concentration, Kj = diag(kf), = diag^kr) and kj,k,. G M"q are non- 
negative forward and reverse kinetic parameters. F, R G are forward and reverse 
stoichiometric matrices. S = — F + R is a consistent stoichiometric matrix defined by the 
existence of at least one strictly positive vector m G M™q, such that S'^ • m = 0. Assuming 
a strictly positive initial concentration xq = x(0) G IR>0' ^'^sn there exists at least one 
non-negative steady state concentration x>q, such that x = 0. 



where f(x) = x{t + r) represents the concentration after an arbitrary small time interval 
r > 0. If it exists, a fixed point x*, such that f(x*) = x*, corresponds to a steady state 
concentration, x*(t) = x*(t + r), or equivalently x = 0. Observe that /(x) is continuous 
as X is given by a continuous function. Let us define the closed, bounded and convex set 



as the domain of f (x). For a strictly positive initial concentration vector, then elementary 
kinetics, continuity of f(x) and non-negative kinetic parameters are sufficient conditions 
to ensure that all subsequent concentrations are non-negative, f(x) > (Theorem 2, [7j). 
Since S G M"^'"' is a stoichiometrically consistent matrix there exists an m G M™q such that 
m-^ -8 = and therefore m-^ • x = 0. It follows that • f(x) = m-^ • xq for all x G 
This together with the non-negativity of f(x) establishes that f(x) G 0,. We have now 
established that f(x) is a continuous mapping from a closed, bounded and convex set into 
itself. By Brouwer's fixed point theorem there exists at least one fixed point f(x*) = x* 
and therefore there exists at least one steady state. □ 

In the proof of existence of steady states for mass conserved elementary kinetics, we 
make use of the following theorem that we state without proof. 

Theorem 2. (Brouwer fixed point theorem). Let Q be a closed, bounded and convex set in 
, and let ^ : M™ — t- be a function that is continuous on $7 and maps into itself. 
Then there exists a point x G such that <I>(x) = x. 

There is a voluminous literature on fixed point theory . Figure IT] illustrates an 



Proof. We define the function f(x) : M™ — > M' 



m 



f (x) = X X 



(5) 



J] = {x > 0, m^ • X = m^ • xo > 0} 
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intuitive appreciation for the rationale behind Brouwer's fixed point theorem by considering 
a one dimensional case, that is, a function mapping of an interval on a line into itself. 

4. Utility of steady state existence theorem 

Theorem [T] is non-constructive in the sense that it does not describe an algorithm for 
computation of steady state concentrations. In the case where one is modeling a system of 
exclusively mass conserved reversible reactions with mass action kinetics, it has long been 
known that a unique steady state concentration can be computed with a single convex 
optimization problem [38j. Such a steady state corresponds to a thermodynamic equilib- 
rium where detailed balance holds. The development of a reliable algorithm to compute 
non-equilibrium steady states for arbitrary large networks is an important open problem. 
In the process of algorithm development, it is essential to know, a priori, if at least one 
steady state exists. Otherwise it becomes impossible to distinguish if a failure to compute 
a steady state is due to a shortcoming of an algorithm's design, or due to to an ill-posed 
problem without a solution in the first instance. 

The key difference between an equilibrium and non-equilibrium steady state is that 
the latter is accompanied by thermodynamic forcing of the system by the environment. In 
chemical reaction networks, time invariant thermodynamic forcing has been mathematically 
represented by clamping a subset of concentrations away from equilibrium or injecting mass 
across the boundary of the model [29]. However, for kinetic modeling of time invariant 
concentrations, any formulation of a forced system must be compatible with the existence 
of at least one steady state concentration vector. It is therefore important to establish, 
if possible, the conditions for existence of at least one steady state for each formulation. 
Some formulations are actually incompatible with the existence of a steady state. 
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4-1- System forcing where a steady state may not exist 

One approach to forcing a system is to represent the exchange of molecules between 
a system and its environment with a set of mass imbalanced source or sink reactions, 
respectively — )■ yl and ^ — )• 0, where A is an arbitrary molecule. Let Se € denote 
the stoichiometry of mass imbalanced exchange reactions. To the author's knowledge, for 
networks with bilinear reactions, no theorem exists that defines the conditions on the data 
{S, Se, kj, k^} such that there still exists at least one non-equilibrium steady state. As 
described in Section j2j| an augmented stoichiometric matrix, S = [ S Sg ]i containing 
mass imbalanced exchange reactions will not be stoichiometrically consistent, so Theorem 
[l] does not apply. 

Another approach to forcing a system, in which all reactions are mass balanced, is to 
attempt to iterate toward a steady state of the forced dynamical system 



where b is a concentration invariant forcing vector in the range of the stoichiometric matrix, 
b G 7^(S). If one chooses a b* G '^(S) such that 



is satisfiable, then this would correspond to forcing in a manner independent of molecule 
concentration. Given x* it is trivial to compute b* but not the other way around. If we as- 
sume that unidirectional reaction rates are as defined in ([T]) and ([2]), then by rearrangement, 
one may express ([3| as 




(6) 




(7) 






12 



where [ F R ] £ Z™"'^". Typically m < n and rank([ F R ]) < so the image of 
[ F R ■ In (x) is not the whole of M^" and therefore the set of all x is a subset of the 
range of the stoichiometric matrix. The set of b* such that ([7]) is satisfiable is only a subset 
of the range of the stoichiometric matrix, so a steady state may not necessarily exist for 
an arbitrary b. Attempting to force a system with ([6]) leaves one with the problem of 
attempting an a priori choice of concentration invariant forcing vector that may not admit 
a steady sate concentration. 

4-2. System forcing where a steady state must exist 

Theorem [T] is constructive in the sense that it leads to a method to force a system in a 
manner that ensures there always exists at least one non-equilibrium steady state concen- 
tration vector. The key point is to recognise that Theorem [T] holds for any choice of non- 
negative kinetic parameters. Additional thermodynamic constraints on kinetic parameters 
([2]) are optional on a per reaction basis. If all reversible reactions have thermodynamically 
feasible kinetic parameters, then the only steady state is thermodynamic equilibrium, but 
if at least one reversible reaction is modeled with thermodynamically infeasible kinetic pa- 
rameters, that violate ([2]), then detailed balance does not hold but there always exists at 
least one non-equilibrium steady state. 

It is not physicochemically realistic to model actual chemical reactions with thermody- 
namically infeasible kinetic parameters for any form of kinetic rate law [8j. However, one 
may include a mass balanced, reversible perpetireaction with thermodynamically infeasible 
kinetic parameters, purely for the modeling purpose of forcing a system away from equilib- 
rium. (The prefix perpeti is from the latin perpes meaning lasting throughout, continuous, 
uninterrupted, continual, perpetual; see www.perseus.tufts.edu.). The augmentation of 
a consistent stoichiometric matrix with a mass balanced perpetireaction still retains the 
stoichiometric consistency of the augmented matrix. Assuming that elementary reaction 
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kinetics is used to model each unidirectional reaction there will still exist a steady state. 
We now illustrate one choice of perpetireaction by considering a biochemical example. 

4-2.1. Mass conserved elementary kinetics of Trypanosoma brucei 

The utility of Theorem [T] can be illustrated by considering a typical kinetic modeling 
scenario, such as the modeling of anaerobic glycolysis in the African trypanosome, Try- 
panosoma brucei, the causative agent of human African trypanosomiasis [3, Ij. Based on a 
phenomenological kinetic model of T. brucei glycolysis [2] the stoichiometry of anaerobic 
glycolysis may be represented in skeleton form by the composite chemical reactions in Fig- 
ure [2j Modeling composite reactions with elementary kinetic rate laws is pseudoelementary 
kinetics, but for the purpose of illustrating the utility of Theorem [T| this distinction is 
superfluous. Starting with extracellular glucose the overall stoichiometry of this pathway 
may be given by the mass balanced composite reaction 

glucose glycerol + pyruvate + . (8) 

This composite reaction may be used as a perpetireaction that, in reverse, connects the 
outputs of anaerobic glycolysis back to the glucose input. This perpetireaction is the 
TrypGlycAner reaction in Figure [2| As perpetireaction kinetic parameters violate Q, no 
equilibrium steady state exists as detailed balance [35] is violated, but there exists at least 
one non-equilibrium steady state for the augmented system. Such a non-equilibrium steady 
state conserves mass in all reactions, and all except the perpetireaction conserve energy. 

Assuming constant temperature and pressure, and uniform spatial concentrations within 
a single compartment, the existence of a single chemical potential for each (compart- 
ment specific) molecule is a necessary and sufficient condition for conservation of energy 
[Ml ESI [30]. The violation of ^ by the pair of forward and reverse elementary kinetic 
parameters, pf and Pr, for the perpetireaction means that there exists no single stan- 
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dard chemical potential for each molecule and hence no single chemical potential for each 
molecule. With reference to the T. hrucei example, one or more of glucose, glycerol, pyru- 
vate or hydrogen ion can not be assigned a unique chemical potential. This is equivalent to 
the statement that the stoichiometrically weighted sum of chemical potential around the 
single stoichiometrically balanced cycle formed by the anaerobic glycolysis pathway and 
the perpetireaction is not zero [U [13] . 

One can also think of the perpetireaction as a chemical reaction that extracts energy, 
but not molecular moieties, from an infinitely large source in the environment. In a non- 
equilibrium steady state, the amount of energy extracted from the environment per unit 
time by the perpetireaction is equal to the entropy production rate of all the other reactions. 
In analogy with electrical networks, at a non-equilibrium steady state, a perpetireaction 
acts like a direct current voltage source. Indeed, in the representation of electrical networks, 
even the most elementary circuit diagram forms a closed cycle with some voltage source 
in the loop to drive electrons around the circuit. In numerically calculated steady states 
the T. hrucei model, if all thermodynamically feasible kinetic parameters are given a unit 
value and Pf > Pr then the net steady state flux of anaerobic glycolysis proceeds in the 
usual direction (Figure [2|. However, if p/ < Pr then anaerobic glycolysis proceeds in the 
reverse direction. In an electrical network analogy, this switch is equivalent to reversing 
the polarity of a direct current voltage source. 

4-2.2. More general examples 

In Section [OTT , we considered an example where the overall system being modeled 



consisted of a single stoichiometrically balanced pathway. More precisely a single extreme 
ray of an augmented stoichiometric matrix 

s=[ s sj 
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where Sp is a column vector with reaction stoichiometry for a single composite perpetire- 
action, with net forward direction right to left in [8} In the general case of an arbitrary 
mass balanced stoichiometric matrix S G M™'", then Sp G M'"''^ can be an arbitrary set of 
k column vectors, each of which could be a linear basis vector for the range of S, accom- 
panied by k perpetireactions. In fact, any b G 7^(8) can be used to create an augmented 
stoichiometric matrix 

S = [ S -b ] 

that is still also consistent. Note however that this is not equivalent to forcing a system 
like |6] as the steady state x*, that we know must exist, satisfies 







-S] •diag( 



k. 



•exp([ F R ]^-ln(x*))-[ b -b ]-diag( 



Pf 

Pr 



exp([ hf hr ]^-ln(x'^)), 



where > are perpetireaction parameters, with bj = max(— b, 0) and b^ = max(b, 0) 

defined in an analogous manner to F and R. Even more general is the consideration 
of continuous kinetic rate laws that guarantee concentration non-negativity. The same 
strategy to augment the system with perpetireactions, yet retain stoichiometric consistency, 
will yield a forced system where there always exists a steady state concentration vector. 

5. Discussion 

In the present work. Theorem [T] gives sufficient conditions for the existence of at least 
one non-negative steady state concentration vector, assuming elementary reaction kinetics 
for a set of mass balanced chemical reactions. All kinetic parameters are required to be 
positive but do not have to satisfy thermodynamic constraints ([2]) on the ratio of forward 
over reverse elementary kinetic parameter. Actual biochemical reactions are modeled with 
reactions that have thermodynamically feasible kinetic parameters. In order to conserve 
mass, yet admit a non-equilibrium steady state, one may augment the set of thermodynam- 
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ically feasible reactions with one or more perpetireactions (perpetual reactions), defined as 
reactions with thermodynamically infeasible kinetic parameters. This gives the flexibility 
to model the non-equilibrium dynamics of a system closed to exchange of mass with the 
environment yet not isolated with respect to the exchange of energy with the environment. 
In a non-equilibrium steady state, the net input of chemical energy is the driving force for 
net flux through the stoichiometrically balanced system. 

Theorem [T] does not preclude that a subset of molecule concentrations are actually zero 
at a steady state. There may exist more than one steady state and no conclusion can be 
drawn as to the stability or otherwise of the steady states that exist. Theorem [T] is non- 
constructive in that it does not provide an algorithm to compute a non-equilibrium steady 
state. However, Theorem [T] makes use of Brouwer's fixed point theorem so, assuming the 
conditions required for Theorem [T] to hold, it may be possible to apply related constructive 
fixed point theorems [15] to design an algorithm that is guaranteed to converge to a non- 
equilibrium steady state. Contributions from fixed point theorists are encouraged and this 



is part of the reason for the detail in Section 2.1 



The sufficient conditions for existence of a non-equilibrium steady state are easily tested 
numerically for arbitrary large chemical networks [32j . As described in an elegant paper by 
Gevorgyan et al. [n] the stoichiometric consistency of a metabolic network can be proved 
or disproved by attempting the linear optimization problem 



minimize e • m (9) 

m 

such that S"^ • m = (10) 

m>0 (11) 

where e denotes a vector of ones and S G ig a stoichiometric matrix. If there exists 



an m satisfying (10) and (11), then S is stoichiometrically consistent, otherwise a suitable 
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solver will provide a certificate of infeasibility indicating that S is inconsistent. Alter- 
natively, if one has rigorously applied mass balancing for each chemical reaction whilst 
reconstructing a network \33\ [3l] , one will be able to assign a positive molecular mass cor- 
responding to each of the molecules in the reconstruction. This strictly positive molecular 



mass vector satisfies (10) and (11), which is sufficient to conclude that the corresponding 
stoichiometric matrix is consistent. 

The only other condition for Theorem [T] to hold is that continuous kinetic rate laws 
be formulated in a manner such that the concentration of any molecule can never be 
negative [7j. In this paper, we have framed Theorem [T] in terms of elementary reaction 
kinetics that satisfy concentration non-negativity if the elementary kinetic parameters are 
non- negative. However, one can envisage a more general version of Theorem [T] as there are 
many other continuous kinetic rate laws that satisfy concentration non-negativity [21J. Any 
continuous kinetic rate law, where the rate of a unidirectional reaction is non-negative and 
zero if and only if any of the concentrations of the molecules consumed in that reaction 
are zero [7j , would form the conditions for a generalised version of Theorem [Tj In any 
case, any phenomenological kinetic rate law can be derived from assumptions that allow 
simplification of a system of elementary chemical reactions |8j. As such, phenomenological 
kinetic modeling has its foundation in mass action kinetics. 

6. Conclusion 

It is 50 years since Wei's Axiom's on the existence of steady states for chemical re- 
action systems |37j and almost 40 years since Horn & Jackson [T7] considered what they 
termed general mass action kinetics. In Horn & Jackson's setting, elementary reaction 
rates are proportional to the abundance of the substrates involved in the reaction, each 
to the power of the absolute value of the corresponding stoichiometric coefficient. How- 
ever, general mass action kinetics considers systems where kinetic parameters need not 
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be thermodynamically feasible, stoichiometric coefficients need not be integers, and mass 
need not be conserved by each reaction. With regard to modehng chemical reaction net- 
works, we agree that consideration of reactions with thermodynamically infeasible kinetic 
parameters does seem profitable for representing the perpetual forcing of a system, purely 
for modeling purposes. However, violation of mass conservation appears unnecessary, as 
a pair of thermodynamically infeasible kinetic parameters are sufficient to force a system 
away from equilibrium, and counterproductive, as the resulting system may not admit a 
non-equilibrium steady state. Horn &: Jackson [17] do realize that the conditions for Wei's 
existence result are unmet when mass is not conserved. We conclude that, rather than 
general mass action kinetics, assuming mass conserved elementary kinetics is sufficient for 
modeling non-equilibrium steady states in arbitrary large biochemical networks, as one is 
then sure that at least one steady state does actually exist. Similar conclusions hold for 
phenomenological kinetic modeling, .e.g., with Michaelis-Menten kinetics, as long as the 
continuos rate laws are such that concentration can never be negative. 

Acknovifledgements 

This work was supported by the U.S. Department of Energy (Offices of Advanced Scien- 
tific Computing Research & Biological and Environmental Research) as part of the Scien- 
tific Discovery Through Advanced Computing program (Grant No. |de-sc0002009). I.T. was 
also supported, in part, by a Marie Curie International Reintegration Grant (No. 249261) 
within the 7th European Community Framework Program. 

References 
References 

[1] Bakker, B. M., Krauth-Siegel, R. L., Clayton, C, Matthews, K., Girolami, M., WesterhofT, H. V., 
Michels, P. A. M., Breitling, R. & Barrett, M. P. 2010. The siUcon trypanosome. Parasitology, 137 
(9), 1333-1341. 



19 



[2] Bakker, B. M., Michels, P. A., Opperdoes, F. R. & Westerhoff, H. V. 1997. Glycolysis in bloodstream 
form trypanosoma brucei can be understood in terms of the kinetics of the glycolytic enzymes. J Biol 

Chcm, 272 (6), 3207-3215. 

[3] Barrett, M. P., Bakker, B. M. & Breitling, R. 2010. Metabolomic systems biology of trypanosomes. 
Parasitology, 137 (9), 1285-1290. 

[4] Beard, D. A., Liang, S. & Qian, H. 2002. Energy balance for analysis of complex metabolic networks. 
Biophys J, 83 (1), 79-86. 

[5] Bernstein, D. S. & Bhat, S. P. 1999. Nonnegativity, reducibility, and semistability of mass action 
kinetics. In Decision and Control, 1999. Proceedings of the 38th IEEE Conference vol. 3, pp. 2206- 
2211, IEEE IEEE. 

[6] Berry, S. R., Rice, S. A. & Ross, J. 2000. Physical Chemistry. 2nd edition, Oxford University Press, 
Oxford. 

[7] Chellaboina, V., Bhat, S., Haddad, W. M. & Bernstein, D. S. 2009. Modeling and analysis of mass- 
action kinetics. Control Systems Magazine, IEEE, 29 (4), 60-78. 

[8] Cook, P. F. & Cleland, W. W. 2007. Enzyme Kinetics and Mechanism. Taylor k. Francis Group, 
London. 

[9] Droste, P., Miebach, S., Niedenfuhr, S., Wiechert, W. & Noh, K. 2011. Visualizing multi-omics data 
in metabolic networks with the software Omix: a case study. Biosystems, 105 (2), 154-161. 

[10] Famili, I. & Palsson, B. 0. 2003. The convex basis of the left null space of the stoichiometric matrix 
leads to the definition of metabolically meaningful pools. Biophys J, 85 (1), 16-26. 

[11] Fell, D. A. & Small, J. R. 1986. Fat synthesis in adipose tissue: an examination of stoichiometric 

constraints. Biochcm J, 238 (3), 781-786. 

[12] Fleming, R. M. T., Maes, C. M., Saunders, M. A., Ye, Y. & Palsson, B. 0. 2012. A variational principle 
for computing nonequilibrium fluxes and potentials in genome-scale biochemical networks. Journal of 
Theoretical Biology, 292, 71-7. 

[13] Fleming, R. M. T., Thiclc, I., Provan, G. & Nashoucr, H. P. 2010. Integrated stoichiometric, thermo- 
dynamic and kinetic modeling of steady state metabolism. J Thoorct Biol, 264, 683-92. 

[14] Gevorgyan, A., Poolman, M. G. & Fell, D. A. 2008. Detection of stoichiometric inconsistencies in 
biomolecular models. Bioinformatics, 24 (19), 2245-2251. 

[15] Granas, A. & Dugundji, J. 2003. Fixed point theory. Springer monographs in mathematics. Springer. 

[16] Heinrich, R. & Schuster, S. 1996. The Regulation of Cellular Systems. Kluwer Academic Publishers, 
Amsterdam. 

[17] Horn, F. & Jackson, R. 1972. General mass action kinetics. Arch Ration Mech Anal, 47 (2), 81-116. 

[18] Jamshidi, N. & Palsson, B. 0. 2008. Formulating genome-scale kinetic models in the post-genome era. 
Mol. Syst. Biol. 4, 171. 

[19] Klipp, E., Liebermeistor, W., Wierling, C, Kowald, A. & Lehrach, H. 2009. Systems biology: a textbook. 
Wiley-Blackwell, Weinheim. 



20 



Lewis, G. 1925. A new principle of equilibrium. Proc Natl Acad Sci USA, 11 (3), 179-183. 

Liebermeister, W., Uhlendorf, J. & Klipp, E. 2010. Modular rate laws for enzymatic reactions: ther- 
modynamics, elasticities and implementation. Bioinformatics, 26 (12), 1528-1534. 

Minty, G. 1960. Monotone networks. Proceedings of the Royal Society of London. Series A, Mathe- 
matical and Physical Sciences, pp. 194 -212. 

Orth, J. D., Thiele, I. & Palsson, B. 0. 2010. What is flux balance analysis? Nat Biotechnol, 28 (3), 
245-248. 

Palsson, B. 0. 2006. Systems Biology: Properties of Reconstructed Networks. Cambridge University 
Press, Cambridge. 

Perelson, A. S. 1976. Remarks on conservation of mass in open chemical reaction systems. J Theor 
Biol, 63 (1), 233-237. 

Planck, M. 1945. Treatise on Thermodynamics. Courier Dover Publications, Chelmsford, MA. 

Price, N. D., Thiele, I. & Palsson, B. 0. 2006. Candidate states of Helicobacter pylori's genome- 
scale metabolic network upon application of "loop law" thermodynamic constraints. Biophys J, 90, 

3919-3928. 

Prigogine, I. & Defay, R. 1954. Chemical thermodynamics. Longmans Green, London. 

Qian, H. & Beard, D. A. 2005. Thermodynamics of stoichiometric biochemical networks in living 
systems far from equilibrium. Biophys Chem, 114 (2-3), 213-220. 

Ross, J. 2008. Thermodynamics and fluctuations far from equilibrium., vol. 80, of Springer Series in 
Chemical Physics. Springer, New York. 

Savinell, J. M. & Palsson, B. 0. 1992. Network analysis of intermediary metabolism using linear 
optimization. I. Development of mathematical formalism. J Theoret Biol, 154 (4), 421-454. 

Thiele, I., Jamsliidi, N., Fleming, R. M. T. & Palsson, B. 0. 2009. Genome-scale reconstruction of E. 
coli's transcriptional and translational machinery: A knowledge-base, its mathematical formulation, 
and its functional characterization. PLoS Comput Biol, 5 (3), el000312. 

Thiele, I. & Palsson, B. 0. 2010. A protocol for generating a high-quality genome-scale metabolic 
reconstruction. Nat Protoc, 5, 93-121. 

Thorleifsson, S. G. & Thiele, I. 2011. rBioNet: A COBRA toolbox extension for reconstructing high- 
quality biochemical networks. Bioinformatics, 27 (14), 2009-2010. 

Tolman, R. C. 1962. The principles of statistical mechanics. Oxford University Press, Oxford, England. 
Watson, M. R. 1986. A discrete model of bacterial metabolism. Comput Appl Biosci, 2 (1), 23-27. 

Wei, J. 1962. Axiomatic treatment of chemical reaction systems. The Journal of Chemical Physics, 36 

(6), 1578-1584. 

White, W., Johnson, S. & Dantzig, G. 1958. Chemical equilibrium in complex mixtures. J Chem Phys, 
28, 751-755. 

Wilhelmy, L. 1850. Ueber das gesetz, nach welchem die einwirkung der sauren auf den rohrzucker 
stattfindet (The law by which the action of acids on cane sugar occurs). Poggendorff's Annalen der 
Physik und Chemie, 81, 413-433. 



21 









^ 1 



Figure 1: Conceptual illustration of the Brouwer fixed point theorem in one dimension. An arbitrary 
continuous function / is represented by the graph mapping the abscissa to the ordinate. On the abscissa, 
the bounded interval between zero and one represents the domain. The domain is closed as it includes 
its endpoints and is convex as every line interval is a convex set. The interval between zero and one on 
the ordinate represents the codomain of /. The image of / is a continuous interval contained within the 
codomain of /, so / is into. It is unnecessary for / to be onto. The diagonal line represents equality of the 
values in the domain and codomain. Tracing a continuous curve from left to right, between the dots, we 
see that it must intersect the diagonal at some point. At the point where the diagonal meets the curve, the 
value of the domain equals the value of the image x* = x' = f{x). When the value passed into the function 
is the same as the value passed out by the function, this value is termed a fixed point. (Figure adapted from 
[http://commons.wikimedia.Org/wiki/File:Fixedpointld.svg) 
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Figure 2: The anaerobic glycolysis pathway of Trypanosoma brucei, part of which is in the cytoplasm 
(c) and part of which is within a membrane-bounded, peroxisome-like organelle termed a glycosome (y). 
The input to this pathway is glucose (glc-D) and the outputs are glycerol (glyc), pyruvate (pyr) and 
hydrogen ion (not shown for clarity). For modeling purposes, one can construct a perpetual reaction, or 
perpetireacUon (TrypGlycAner) , from extracellular output metabolites to extracellular glucose input to 
form, with the anaerobic glycolysis pathway, a stoichiometrically balanced cycle. With appropriate choice 
of kinetic parameters, a non-equilibrium steady state concentration corresponds to net flux in the directions 
indicated by the arrows. Dashed arrows indicate the involvement of cofactors. (Illustration created with 
Omix j9j.) 
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